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Atomistically  Informed  Continuum  Model  of  Polymer-Based  Nanoeomposites 

Catalin  R.  Picu,  Alireza  Sarvestani  and  Murat  S.  Ozmusul 
Department  of  Mechanical,  Aerospace  and  Nuclear  Engineering 
Rensselaer  Polytechnic  Institute,  Troy,  NY  12180 

ABSTRACT 

A model  polymeric  material  filled  with  spherical  nanoparticles  is  considered  in  this  work. 
Monte  Carlo  simulations  are  performed  to  determine  the  polymer  chain  conformations  in  the 
vicinity  of  the  curved  interface  with  the  filler.  Several  discrete  models  of  increasing  complexity 
are  considered:  the  athermal  system  with  excluded  volume  interactions  only,  the  system  in  which 
entropic  and  energetic  interactions  take  place  while  the  filler  is  a purely  repulsive  sphere,  and  the 
system  in  which  both  filler-polymer  and  polymer-polymer  energetic  interactions  are  accounted 
for.  The  total  density,  chain  end  density,  chain  segment  preferential  orientation  and  chain  size 
and  shape  variation  with  the  distance  from  the  filler  wall  are  determined.  The  structure  is  graded, 
with  the  thickness  of  the  transition  region  being  dependent  on  the  property  and  scale  considered. 
Hence,  the  polymer  in  the  vicinity  of  the  filler  is  represented  in  the  continuum  sense  by  a 
material  with  graded  properties  whose  elasticity  is  determined  based  on  the  local  structure. 
Homogenization  theory  is  the  used  to  obtain  the  overall  composite  moduli.  The  filler  size  effect 
on  the  composite  elasticity  is  evaluated. 

INTRODUCTION 

Polymer-based  nanocomposites  form  a new  class  of  reinforced  polymers  in  which  the  fillers 
have  at  least  one  dimension  less  than  100  nm.  These  materials  have  macroscopic  properties 
significantly  different  from  those  of  conventional  composites  of  same  filler  volume  fraction 
[e.g.1,2].  For  example,  the  dispersion  of  only  2%  volume  fraction  of  nanoparticles  in 
thermoplastics  almost  doubles  the  yield  stress  and  increases  the  Young’s  modulus. 

The  mechanisms  responsible  for  the  enhanced  properties  are  not  fully  understood  and 
controlled.  It  appears  that  interesting  properties  are  obtained  when  the  filler  size  and/or  the  filler 
wall-to-wall  distance  become  comparable  with  the  chain  radius  of  gyration.  Two  principal  ideas 
have  been  promoted  in  this  connection:  the  novel  properties  are  either  due  to  chains  connecting 
several  fillers  (the  double  network  theory),  or  to  the  special  distribution  and  density  of 
entanglements  in  the  confined  polymer  matrix.  In  all  situations,  it  is  generally  accepted  that  the 
structure  of  polymeric  chains  confined  between  fillers  and  their  binding  to  the  wall  are  the  major 
elements  controlling  the  properties  of  the  composite. 

This  work  is  part  of  an  ongoing  effort  to  relate  the  polymer  structure  to  the  macroscopic 
nanocomposite  properties.  We  consider  a model  polymer  filled  with  spherical  impenetrable 
particles.  The  chain  structure  in  the  vicinity  of  fillers  is  determined  by  Monte  Carlo  simulations. 
The  elasticity  of  the  matrix  is  determined  based  on  the  local  structure.  The  composite  is 
homogenized  for  its  overall  elastic  properties.  We  first  review  the  main  features  of  the  polymer 
structure  at  the  wall,  then  we  present  the  procedure  by  which  the  continuum  model  is  calibrated 
based  on  the  discrete  structure,  and  finally,  we  discuss  the  filler  size  effect  on  composite  moduli. 
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Polymer  structure  in  the  vicinity  of  fillers 


Lattice  Monte  Carlo  simulations  are  performed  to  determine  the  polymer  structure  in  the 
vicinity  of  fillers.  A simulation  cell  containing  a single  fdler  is  used  (with  periodic  boundary 
conditions).  The  polymers  are  represented  by  a “pearl-necklace”  model.  The  beads  are  free  to 
move  on  a lattice  surrounding  the  filler,  the  chain  connectivity  and  the  excluded  volume 
conditions  being  enforced  at  all  times.  The  model  and  the  simulation  procedure  are  described  in 
detail  in  [3].  The  parameters  of  the  model  are  the  chain  size  (N),  the  filler  radius  (R)  and  the 
filler  volume  fraction.  Various  types  of  interactions  are  considered  in  separate  models:  excluded 
volume  only,  excluded  volume  and  energetic  interaction  in  bulk  polymer,  excluded  volume  and 
energetic  interaction  between  polymers  and  the  filler.  The  results  reported  here  arc  for  the  system 
with  R = 8,  N = 100  (Rg  = 4.7)  and  filler-to-filler  distance  L = 32  (wall-to-wall  distance  of  16, 
i.e.  twice  the  chain  size,  4Rg).  All  dimensions  are  in  multiples  of  the  bond  length  (equal  to  the 
unit  of  length  used  in  the  Lennard-Jones  potential  representing  non-bonded  interactions). 

The  total  bead  number  density  normalized  by  the  bulk  density  is  shown  in  Fig.  la  for  the  case 
in  which  no  filler-polymer  interactions  are  imposed.  The  curve  labeled  “athermal”  corresponds 
to  the  purely  entropic  system  in  which  only  excluded  volume  interactions  are  allowed.  This 
corresponds  to  the  “high  temperature”  situation  in  which  the  entropic  component  of  the  free 
energy  dominates.  The  curve  labeled  “energetic”  corresponds  to  the  system  in  which  cohesive 
interactions  in  the  bulk  polymer  are  considered.  An  ample  depiction  layer  forms  next  to  the  filler, 
and  effect  also  observed  in  lattice  MC  simulations  of  structure  next  to  flat  interfaces.  We  note 
that  continuous  space  discrete  models  predict  strong  density  oscillations  next  to  the  wall  due  to 
packing  and  no  depleted  layer.  The  effect  seen  here  is  pronounced  in  the  energetic  case  (low 
temperatures),  while  in  the  athermal  system  it  is  absent,  the  bead  density  being  constant. 

An  entropic  driving  force  exists  for  chain  end  segregation  in  the  interface.  The  chain  end 
number  density  is  shown  in  Fig  la  for  the  energetic  system.  The  end  segregation  has  no  effect  on 
the  elastic  moduli  of  the  matrix.  The  effect  of  the  various  model  parameters  is  discussed  in  [3). 

The  chain  number  density  is  evaluated  based  on  the  position  of  the  chain  centers  of  mass 
with  respect  to  the  wall  (Fig.  lb).  The  volume  surrounding  the  filler  is  divided  in  spherical  bins 
(concentric  layers)  of  thickness  equal  to  2.  The  data  points  represent  the  density  in  each  layer 
normalized  with  the  bulk  chain  number  density.  The  chain  density  is  lower  at  the  wall  than  in 
the  bulk,  an  effect  due  to  the  excluded  volume  of  the  filler.  The  depletion  is  also  driven  by  the 


distance  from  filler  surface  distance  from  filler  surface 

Figure  1.  Bead  (a)  and  chain  (b)  number  density  profile  next  to  the  spherical  filler. 
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entropic  force  that  retracts  the  chains  from  the  wall  (spherical  or  planar),  and  by  the  densification 
effect  of  the  cohesive  interactions  in  the  bulk  polymer  (see  the  difference  between  the 
“energetic”  and  “athermal”  curves),  A small  number  of  chains  that  wrap  around  the  filler  have 
their  centers  of  mass  within  the  particle  (negative  distance  from  filler  surface). 

The  chain  size  was  determined  by  computing  the  gyration  tensor  for  chains  having  their 
centers  of  mass  at  various  distances  from  the  filler  surface.  In  the  bulk,  chains  have  an  ellipsoidal 
shape,  the  semi-axes  of  the  ellipsoid  being  given  by  the  eigenvalues  of  the  gyration  tensor.  The 
results  obtained  for  the  present  system  are  shown  in  Fig.  2.  The  three  semi-axes  of  magnitude  in 
the  bulk  X]  = 3.89,  = 1 .84  and  1.14,  are  shown  in  Fig.  2a.  Interestingly,  chains  having 

their  centers  of  mass  in  the  first  bin  next  to  the  wall  have  dimensions  similar  to  those  in  the  bulk. 
The  few  chains  having  their  centers  of  mass  within  the  filler  are  swollen. 

The  orientation  of  the  long  semi-axis  of  the  ellipsoidal  chains  is  random  in  the  bulk.  As  the 
chain  center  of  mass  approaches  the  wall,  the  chains  tend  to  align  with  their  large  semi-axis  in 
the  direction  tangential  to  the  interface  (a  “docking”  transition  previously  observed  next  to  a flat 
wall).  Fig.  2b  shows  the  second  moment  of  the  orientation  distribution  of  the  long  semi-axis 
(labeled  “chain”).  Perfect  orientation  in  the  tangential  direction  corresponds  to  a value  of  -0.5  of 
this  measure.  Preferential  orientation  is  observed  on  all  scales,  from  the  bond  to  the  chain  scale. 
The  orientation  of  bonds  and  chain  end-to-end  vectors  is  weaker  than  that  of  the  coil  (Fig.  2b). 

Higher  confinement  induced  by  a further  reduction  of  the  wall-to-wall  distance  leads  to  a 
reduction  of  chain  dimensions,  with  chains  “docking”  to  the  wall  in  all  cases  [3].  This  result  is  in 
agreement  with  data  published  by  Vacatello  [4],  but  disagree  with  the  model  advanced  by  Mark 
and  collaborators  [5],  Our  data  are  discussed  against  both  sets  of  results  in  [3], 

Elastic  moduli  of  the  system  of  chains 


The  elasticity  of  the  system  of  chains  is  derived  here  based  on  the  free  energy  of  the  system. 
The  total  free  energy,  W,  contains  three  energetic  and  an  entropic  term: IP  =Unh  + Uh  +U  f -TS 
The  energetic  terms  are  due  to  non-bonded  interactions,  between  beads  not  directly  connected 
along  a chain,  to  bonded  interactions,  between  directly  connected  beads,  and  to  the  energetic 
interaction  between  the  filler  and  the  matrix.  The  stiffness  of  the  system  of  chains  in  the  model 


distance  from  filler  surface  distance  from  filler  surface 

Figure  2.  (a)  Semi-axes  of  the  ellipsoidal  chains  as  a function  of  the  position  of  their  centers 
of  mass  with  respect  to  the  wall,  (b)  A measure  of  the  preferential  orientation  of  chains,  bonds 
and  chain  end-to-end  vector  in  the  direction  tangential  to  the  spherical  filler. 
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results  as  the  second  derivative  of  the  free  energy  density  (w  = W/V)  with  respect  to  the 
deformation  gradient  tensor  F,  CjJkI  = 92w/3Fi,dFkr.  Hence,  the  stiffness  has  4 components 
corresponding  to  the  respective  terms  in  the  free  energy. 

The  total  energy  due  to  non-bonded  interactions,  Unb , reads 

U"b  = Ju"bdVm  = J~  |n (rRA , 0RA , (pRA  )<J»l,b(rRA  )d VrdVm  , (1) 

vm  Vr 

where  unb  represents  the  non-bonded  energy  density,  and  V„,  is  the  volume  of  the  polymer  matrix 
within  the  simulation  cell.  The  non-bonded  energy  density  may  be  expressed  as  the  energy  due  to 
non-bonded  interactions  (cohesive  energy)  of  a representative  atom,  RA,  divided  by  the  atomic 
volume  Q.  In  the  integral  expression  of  this  term,  <J)nh  is  the  interatomic  potential  and  n"b  stands 
for  the  density  of  non-bonded  neighbors  of  the  RA.  The  energy  Unb  is  estimated  here  within  the 
local  approximation.  The  inner  integral  in  eqn.  (1)  is  evaluated  by  neglecting  the  dependence  of 
n"b  on  (0RA , <pRA ) , i.e.  considering  that  an  atom  located  in  bin  p has  a neighborhood  similar  to 
that  of  an  RA  located  in  a homogeneous  bulk  of  density  p(rp).  where  rp  denotes  the  radial 
position  of  bin  p with  respect  to  the  coordinate  system  centered  at  the  center  of  the  filler.  The 
function  p(rp)  is  shown  in  Fig.  la.  For  an  atom  located  in  a homogeneous  material  of  uniform 
density  p,  the  function  n"b  may  be  expressed  as  n”b  = p g(rRA),  where  g(rRA>  is  the  pair 
distribution  function  of  non-bonded  neighbors  about  the  RA.  Then,  the  energy  density  due  to 
non-bonded  interactions  in  bin  p may  be  written 

un„.p  = 1 A-  fpb«,kg(rRA)<t,nb(rRA)dVr  = (2) 

& Phulk  V,  P bulk 

Hence,  the  local  stiffness  in  bin  p,  associated  with  non-bonded  interactions,  results  from  the 
stiffness  of  the  bulk  by  scaling  with  the  ratio  of  the  local  and  bulk  densities 
C,"b;p  = (pr  /pblllk  )C|jbb“'k  Since  the  density  varies  with  the  distance  from  the  wall,  this 
relationship  suggests  that  the  matrix  may  be  represented  in  the  continuum  sense  by  a material 
with  graded  elastic  constants. 

The  contribution  to  stiffness  of  bonded  interactions  may  be  evaluated  by  a similar  procedure. 
In  the  model  considered  here,  the  bonds  arc  assumed  to  be  rigid  and  hence,  Ub  = 0. 

No  energetic  interactions  between  the  filler  and  polymers  are  considered.  Hence.  Ub  = 0.  The 
more  realistic  case  in  which  the  polymers  arc  attracted  to  the  filler,  while  cohesive  interactions 
are  enabled  in  the  bulk,  is  discussed  in  [6]. 

The  entropic  contribution  to  stiffness  may  be  evaluated  based  on  the  classical  formulation  of 
stress  production  in  polymeric  systems  and  the  entropic  spring  idea.  The  entropy  of  a chain  is 
related  to  the  probability  distribution  of  the  end-to-end  vector.  R'e.  For  Gaussian  chains  (an 
assumption  shown  to  be  valid  in  the  system  discussed  here),  the  entropic  contribution  to  the  free 

2 ■> 

energy  per  chain  reads  W = (3/2)kTRcc  /(C,b  N) . A deformation  described  by  the  deformation 
gradient  tensor  F changes  the  end-to-vector  from  R“  in  the  undeformed  configuration  into 
R“  = F • R“  . Then,  the  contribution  of  a chain  to  the  stiffness  tensor  becomes 
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(3) 
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where  pc  is  the  chain  number  density  (given  in  Fig.  lb),  and  6;s,  is  the  Kroneeker  delta.  The  total 
entropic  contribution  to  the  stiffness  results  by  integrating  over  the  whole  chain  orientation 
distribution.  In  the  bulk,  the  chains  are  randomly  oriented  in  space  and  the  entropic  contribution 
is  given  by  Cf£[  = p'kTS^S^ . In  the  neighborhood  of  the  filler,  the  end-to-end  chain  vector  is 
preferentially  oriented  in  the  direction  tangential  to  the  filler  (Fig.  2b)  which  leads  to  an 
anisotropic  C“k" . As  discussed  in  [6],  the  anisotropy  is  rather  weak  and  hence  is  neglected  in  this 
analysis.  It  results  that  the  entropic  contribution  to  stiffness  in  bin  p is  proportional  to  the 
entropic  contribution  in  the  bulk,  the  proportionality  factor  being  the  ratio  of  the  chain  number 
density  in  bin  p to  that  in  the  bulk, 

/-ientr,p  / c / c \/-<entr,bulk 

'“'iJkL  \Pp  / Pbulk/'-'UkL 

Hence,  the  stiffness  in  bin  p normalized  by  the  bulk  stiffness  reads 

Op 

Hjk!. 

pbulk  VKp/Kbulk 


= (Pp/Pbuik)Qi  +(p;/pU)5it6JLQ2,  where  Q,  =C*** /Cjj  and  Q2  =C"/C“. 


The  variation  of  the  two  density  ratios  with  the  distance  from  the  wall  is  shown  in  Figs,  la  and  b, 
while  Qj  and  Q2  are  constants,  independent  of  the  distance  to  the  wall.  The  two  terms  Qi  and  Q2 
represent  the  fraction  of  the  total  stiffness  of  the  bulk  due  to  energetic  and  entropic  interactions, 
respectively.  Their  magnitude  depends  on  temperature  and  bulk  density.  Here,  we  take  these 
quantities  as  parameters  (note  that  Qi  + Qi  = 1).  At  low  temperatures,  the  entropic  component  is 
small  and  the  first  term  dominates.  As  the  temperature  increases,  the  second  term  becomes 
important.  The  stiffness  gradient  in  the  neighborhood  of  the  filler  corresponding  to  3 pairs 
(Qi.Ch)  is  shown  in  Fig.  3a. 

Once  the  stiffness  tensor  is  defined  in  each  layer,  homogenization  theory  may  be  used  to 
evaluate  the  moduli  of  the  composite.  For  this  purpose,  the  multi-inclusion  model  proposed  by 
Nemat-Nasser  and  Hori  [7]  was  used.  The  expression  for  the  overall  composite  stiffness  tensor  is 
given  by  C = Cb,m  [I  + (S  - 1)  0>R  ] (I  + S Or  )‘ 1 , where  Cbulk  is  the  stiffness  tensor  of  the  bulk, 
S(fl;Cbulk)  is  Eshelby’s  tensor  for  a spherical  inclusion  <T2  embedded  in  a solid  of  stiffness  C1 


(bulk 


and  Or  = £(fp -fp-|)'F!’(£2;CMk,Cp),  yp  = -[S(f2;CMk)  + (Cp -Cbu"T 


C”]'1  [7], 
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Filler  size  effect 


Decreasing  the  size  of  the  filler  (R)  at  constant  volume  fraction  requires  a proportional 
reduction  of  the  cell  size  and  of  the  wall-to-wall  distance  between  fillers.  In  this  process,  the 
thickness  of  the  graded  layer,  which  is  essentially  a function  of  the  chain  length  N,  does  not 
scale.  Hence,  the  volume  fraction  occupied  by  the  graded  material  increases  and  a size  effect  in 
the  overall  composite  moduli  is  evidenced.  The  bulk  and  shear  moduli  of  the  composite  are 
evaluated  for  the  system  with  6.5%  filler  volume  fraction.  The  result  for  the  shear  modulus  is 
shown  in  Fig.  3b  for  three  (Qi,  Q>)  pairs.  The  composite  modulus  is  normalized  by  the  modulus 
of  the  corresponding  conventional  composite,  with  no  graded  region,  and  having  the  same  filler 
volume  fraction.  In  the  limit  of  large  fillers,  the  modulus  tends  to  that  of  the  conventional 
composite  (to  1 in  Fig.  3b).  As  the  minimum  wall-to-wall  distance  between  fillers  becomes 
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distance  from  filler  surface  filler  radius 


Figure  3.  (a)  Variation  of  matrix  shear  modulus  with  the  distance  from  the  wall.  The  variation  is 
mainly  due  to  reduced  bead  and  chain  number  density  at  the  wall,  (b)  Scaling  of  the  composite 
shear  modulus  with  the  size  of  the  filler  (R)  at  constant  filler  volume  fraction  (6.5%). 

smaller  than  about  5Rg,  the  modulus  departs  from  that  of  the  conventional  composite  exhibiting  a 
sensible  size  effect.  The  minimum  filler  separation  considered  was  about  2Rg  (which 
corresponds  to  R = 4)  in  order  for  the  graded  regions  of  neighboring  fillers  not  to  overlap.  The 
shear  modulus  decreases  as  the  filler  size  is  reduced.  The  bulk  modulus  has  a similar  variation. 

CONCLUSIONS 

The  structure  of  the  polymer  matrix  in  the  close  neighborhood  of  a spherical  nanofiller  in 
polymer-based  nanocomposites  is  studied.  It  is  shown  that  the  chains  having  their  center  of  mass 
close  to  the  filler,  rotate  without  deforming  with  their  large  semiaxis  in  the  direction  tangential  to 
the  filler.  The  chain  number  density  is  reduced  at  the  wall.  The  stiffness  tensor  of  this  graded 
material  was  derived  based  on  the  free  energy  of  the  system  which,  in  turn,  was  evaluated  based 
on  the  atomic  and  molecular  structure.  The  composite  was  represented  by  a multi-phase 
continuum  model  and  homogenized  using  the  homogenization  theory.  The  effect  on  the  overall 
composite  moduli  of  reducing  the  filler  size  at  constant  filler  volume  fraction  was  evaluated. 
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